#################################################################################
# Replication file for: "Accounting for Noncompliance in Survey Experiments"    #
#                                                                               #
# Jeffrey J. Harden                                                             #
# University of Notre Dame                                                      #
# jeff.harden@nd.edu                                                            #
#                                                                               #
# Anand E. Sokhey                                                               #
# University of Colorado Boulder                                                #
# anand.sokhey@colorado.edu                                                     #
#                                                                               #
# Katherine L. Runge                                                            #
# University of Colorado Boulder                                                #
# katherine.runge@colorado.edu                                                  # 
#                                                                               #
# Additional replications                                                       #
# Last update: August 2, 2018                                                   #
#################################################################################
### Packages ###
library(foreign)
library(AER)
library(ggplot2)

### Replications ###
## Corbacho et al. (2016, AJPS) [Meta Analysis ID 20] ##
# Note: This paper reports the ITT and CACE, so we only
# need to run the original replication files.
source("corbacho-1.R")
corba.itt1 <- ITT
corba.itt.se1 <- ITT.se

source("corbacho-2.R")
corba.itt2 <- ITT
corba.itt.se2 <- ITT.se

source("corbacho-3.R")
corba.cace1 <- LATE.corruption[1]
corba.cace.se1 <- LATE.corruption[2]

source("corbacho-4.R")
corba.cace2 <- LATE.corruption[1]
corba.cace.se2 <- LATE.corruption[2]
  
# Compile results #
corba.res <- data.frame(article = "Corbacho et al. (2016)",
                        estimate = c("corruption no covariates", "corruption with covariates"),
                        itt = c(corba.itt1, corba.itt2),
                        cace = c(corba.cace1, corba.cace2),
                        est.ratio = c(corba.cace1, corba.cace2)/c(corba.itt1, corba.itt2),
                        itt.se = c(corba.itt.se1, corba.itt.se2),
                        cace.se = c(corba.cace.se1, corba.cace.se2),
                        se.ratio = c(corba.cace.se1, corba.cace.se2)/c(corba.itt.se1, corba.itt.se2))
rownames(corba.res) <- NULL

## Grose et al. (2015, AJPS) [Meta Analysis ID 25] ##
# Note: The manipulation check question does not appear to be in the replication data.
# This replication uses latency, with a cutoff of 100 seconds (the mean reading time),
# to measure compliance.
# See the Stata files table4-itt.do and table4-cace.do for replication code.

gmvh1 <- read.csv("table4-itt.csv")
gmvh2 <- read.csv("table4-cace.csv")

# Compile results #
gmvh.res <- data.frame(article = "Grose et al. (2015)",
                        estimate = paste(gmvh1$DV, gmvh1$Treatment, sep = "."),
                        itt = gmvh1$Mean,
                        cace = gmvh2$Mean,
                        est.ratio = gmvh2$Mean/gmvh1$Mean,
                        itt.se = gmvh1$Tse,
                        cace.se = gmvh2$Tse,
                        se.ratio = gmvh2$Tse/gmvh1$Tse)
rownames(gmvh.res) <- NULL
gmvh.res$est.ratio[12] <- 1 # Possible problematic CACE estimate

## Ahler (2014, JOP) [Meta Analysis ID 43] ##
# Data #
ahler <- read.dta("study_2.dta")

# Compliance variable (defined in online appendix, page 34) ##
ahler$comply.distort_t <- ahler$comply.tell_t <- 0
ahler[ahler$manipcheck %in% c(0, 2) & ahler$distort_t == 1, ]$comply.distort_t <- 1
ahler[ahler$manipcheck == 1 & ahler$tell_t == 1, ]$comply.tell_t <- 1

# Estimate ITT and CACE #
ahler.itt <- lm(moderateness ~ distort_t + tell_t, data = ahler)
ahler.cace <- ivreg(moderateness ~ comply.distort_t + comply.tell_t | distort_t + tell_t, data = ahler) 

# Compile results #
ahler.res <- data.frame(article = "Ahler (2014)",
                        estimate = c("distort condition", "tell condition"),
                        itt = coef(ahler.itt)[2:3],
                        cace = coef(ahler.cace)[2:3],
                        est.ratio = abs(coef(ahler.cace)[2:3])/abs(coef(ahler.itt)[2:3]),
                        itt.se = sqrt(diag(vcov(ahler.itt))[2:3]),
                        cace.se = sqrt(diag(vcov(ahler.cace))[2:3]),
                        se.ratio = sqrt(diag(vcov(ahler.cace))[2:3])/sqrt(diag(vcov(ahler.itt))[2:3]))
rownames(ahler.res) <- NULL

## Horowitz and Levendusky (2011, JOP) [Meta Analysis ID 52] ##
# Data #
hl <- read.dta("conscription_data.dta", convert.factors = FALSE)

# Outcome variables
hl$send.troops <- ifelse(hl$q3b1_2 > 0, (-1*hl$q3b1_2 + 5), NA)
hl$binary.send.troops <- ifelse(hl$send.troops > 2, 1, 0)

# Covariates
hl$hawk <- ifelse(hl$q1b > 0, hl$q1b, NA)
hl$military <- ifelse(hl$pppa0045 > 0, ifelse(hl$pppa0045 < 4, 1, 0), NA)
hl$female <- ifelse(hl$ppgender == 2, 1, 0)
hl$age <- hl$ppage/10
hl$income <- hl$ppincimp
hl$college <- ifelse(hl$ppeducat == 4, 1, 0)
hl$pid <- (-1*hl$xparty7) + 8 

# Experimental conditions 
hl$saw.draft.info <- ifelse(hl$xtessp2 == 1 | hl$xtessp2 == 3, 1, 0)
hl$high.casualties <- ifelse(hl$xtessp2 == 1 | hl$xtessp2== 2, 1, 0)

# Compliance variable (defined on page 529) ##
hl$draft.manip.check <- ifelse(hl$q4b > 0, hl$q4b, NA)
hl$comply <- 0
hl$sdi.comply <- ifelse(hl$saw.draft.info == 1 & 
                        hl$draft.manip.check >= 3, 1, 0)

# Estimate ITT and CACE #
hl.itt <- lm(binary.send.troops ~ high.casualties*saw.draft.info + 
                                  military + female + age + income + 
                                  college + pid + hawk, data = hl)
hl.cace <- ivreg(binary.send.troops ~ high.casualties*sdi.comply + 
                                      military + female + age + income + 
                                      college + pid + hawk | 
                                      high.casualties*saw.draft.info + 
                                      military + female + age + income + 
                                      college + pid + hawk, data = hl)

# Compile results #
hl.res <- data.frame(article = "Horowitz and Levendusky (2011)",
                     estimate = c("high casualties condition",
                                  "saw draft information condition",
                                  "high casualties x saw draft condition"),
                     itt = coef(hl.itt)[c(2, 3, 11)],
                     cace = coef(hl.cace)[c(2, 3, 11)],
                     est.ratio = abs(coef(hl.cace)[c(2, 3, 11)])/abs(coef(hl.itt)[c(2, 3, 11)]),
                     itt.se = sqrt(diag(vcov(hl.itt))[c(2, 3, 11)]),
                     cace.se = sqrt(diag(vcov(hl.cace))[c(2, 3, 11)]),
                     se.ratio = sqrt(diag(vcov(hl.cace))[c(2, 3, 11)])/sqrt(diag(vcov(hl.itt))[c(2, 3, 11)]))
rownames(hl.res) <- NULL

## Stephens-Dougan (2016, JOP) [Meta Analysis ID 64] ##
# Data #
d.sd <- read.dta("sd.dta")

# Compliance variable (defined on page 691) ## 
d.sd$comply <- ifelse(d.sd$manipchk_environ == "no" &
                      d.sd$manipchk_black == "no" &
                      d.sd$manipchk_afghanistan == "no" &
                      d.sd$manipchk_glasses == "no" &
                      d.sd$manipchk_na == "yes", 1, 0)

d.sd$wdw.comply <- ifelse(d.sd$whitedem_whites == 1 & d.sd$comply == 1, 1, 0)       
d.sd$wdd.comply <- ifelse(d.sd$whitedem_diverse == 1 & d.sd$comply == 1, 1, 0)  
d.sd$wdb.comply <- ifelse(d.sd$whitedem_blacks == 1 & d.sd$comply == 1, 1, 0)  
d.sd$wrw.comply <- ifelse(d.sd$whiterepub_whites == 1 & d.sd$comply == 1, 1, 0)  
d.sd$wrd.comply <- ifelse(d.sd$whiterepub_diverse == 1 & d.sd$comply == 1, 1, 0)  
d.sd$wrb.comply <- ifelse(d.sd$whiterepub_blacks == 1 & d.sd$comply == 1, 1, 0)  

# Estimate ITT and CACE #
sd.outcomes <- c("votegregdavis", "fair_gregdavis", "gregdavis_redcrim",
                 "davis_affirmact", "blacksoverwhites")
sd.est.itt <- sd.est.cace <- sd.se.itt <- sd.se.cace <-  matrix(NA, nrow = 5, ncol = 5)
rownames(sd.est.itt) <- rownames(sd.est.cace) <- c("whitedem_diverse", "whitedem_blacks",
                                                   "whiterepub_whites", "whiterepub_diverse",
                                                   "whiterepub_blacks")
colnames(sd.est.itt) <- colnames(sd.est.cace) <- colnames(sd.se.itt) <- colnames(sd.se.cace) <- sd.outcomes

for(i in 1:length(sd.outcomes)){
sd.itt <- lm(as.formula(paste(sd.outcomes[i], "~ whitedem_diverse + whitedem_blacks + 
                                              whiterepub_whites + whiterepub_diverse + 
                                              whiterepub_blacks + gender + educ +
                                              income + pid7 + south")),
                                              weights = weight, data = d.sd)
sd.cace <- ivreg(as.formula(paste(sd.outcomes[i], "~ wdd.comply + wdb.comply + wrw.comply +
                                                     wrd.comply + wrb.comply + gender +
                                                     educ + income + pid7 + south |
                                                     whitedem_diverse + whitedem_blacks +
                                                     whiterepub_whites + whiterepub_diverse +
                                                     whiterepub_blacks + gender + educ +
                                                     income + pid7 + south")),
                                                     weights = weight, data = d.sd)
sd.est.itt[i, ] <- coef(sd.itt)[2:6]
sd.est.cace[i, ] <- coef(sd.cace)[2:6]
sd.se.itt[i, ] <- sqrt(diag(vcov(sd.itt))[2:6])
sd.se.cace[i, ] <- sqrt(diag(vcov(sd.cace))[2:6])
}

# Compile results #
nms <- NULL
for(i in 1:nrow(sd.est.itt)){
  for(j in 1:ncol(sd.est.itt)){
    nms <- c(nms, paste(rownames(sd.est.itt)[i], sd.outcomes[j], sep = "."))
  }
}

sd.res <- data.frame(article = "Stephens-Dougan (2016)",
                     estimate = nms,
                     itt = c(t(sd.est.itt)),
                     cace = c(t(sd.est.cace)),
                     est.ratio = abs(c(t(sd.est.cace)))/abs(c(t(sd.est.itt))),
                     itt.se = c(t(sd.se.itt)),
                     cace.se = c(t(sd.se.cace)),
                     se.ratio = abs(c(t(sd.se.cace)))/abs(c(t(sd.se.itt))))
rownames(sd.res) <- NULL

## Harden (2016) ##
# Note: This replication uses a latency of 35 seconds as the compliance cutoff.
load("harden-replication.RData")

harden.res <- data.frame(article = "Harden (2016)",
                         estimate = "pork barrel",
                         itt = coef(itt)[2],
                         cace = coef(cace35)[2],
                         est.ratio = abs(coef(cace35)[2])/abs(coef(itt)[2]),
                         itt.se = sqrt(diag(vcov(itt))[2]),
                         cace.se = sqrt(diag(vcov(cace35))[2]),
                         se.ratio = abs(sqrt(diag(vcov(cace35))[2]))/abs(sqrt(diag(vcov(itt))[2])))
rownames(harden.res) <- NULL

### Summarize Results ###
final.res <- rbind(corba.res, gmvh.res, ahler.res, hl.res, sd.res, harden.res)

theme_set(theme_gray(base_size = 18))
x.la <- paste(seq(0, 250, 25), "%", sep = "")

## Estimates ##
# Figure A3(a) #
pdf("figureA3a.pdf")

ggplot(final.res, aes(x = est.ratio)) +
  geom_histogram(bins = 35) +
  scale_x_continuous(breaks = seq(1, 3.5, .25), labels = x.la) +
  scale_y_continuous(breaks = 0:15) +
  xlab("% increase in magnitude, ITT to CACE") +
  ylab("Count") +
  annotate("text", x = 2.6, y = 4, size = 5,
           label = "Median: 28%    \n IQR: (5%, 71%)")

dev.off()

## SEs ##
# Figure A3(b) #
pdf("figureA3b.pdf")

ggplot(final.res, aes(x = se.ratio)) +
  geom_histogram(bins = 35) +
  scale_x_continuous(breaks = seq(1, 3.5, .25), labels = x.la) +
  scale_y_continuous(breaks = 0:15) +
  xlab("% increase in magnitude, ITT to CACE") +
  ylab("Count") +
  annotate("text", x = 2.6, y = 4, size = 5,
           label = "Median: 59%      \n IQR: (47%, 84%)")

dev.off()

### Save Workspace ###
save.image("additional-replications.RData")







